% Run Counterfactuals
%
%

clear all
clc

warning('Off','all') ;

path(pathdef) ;

addpath('./routines')
addpath('./routines/solmat')

%% Load posterior chain

load mhall_01-Apr-2021

SET.EST.burn_in = 1/2;

params_x_ = [] ;
 
for runs=2:maxproc
    params_x_ = [params_x_ ; params_x(round(length(params_x)*SET.EST.burn_in):end,:,runs)] ; %#ok<AGROW>
end

for ii=1:length(params_x_(1,:))
    [a,b]=ksdensity(params_x_(:,ii),'NumPoints',200);
    est_params(ii) = b(a==max(a)) ;
end

save est_params est_params

%% Run Dynare

clear
clc

setup_model

%% Load and setup demographic parameters

trend_dates = 1940:0.25:2070 ; 

%convert to quarterly
year_vec  = 1940:2142 ;
datevec_q = 1940:0.25:2142 ;

start_date = 1986 ; SET.start_date=start_date;
end_date   = 2019.75; SET.end_date=end_date;

trend_start = 1960 ;
trend_end   = 2060 ;

trend_dates = trend_dates(find(trend_dates==trend_start):find(trend_dates==trend_end)) ;
SET.trend_dates = trend_dates ;

SET.usecore = 0 ; 
SET.unant_demographics = 0 ;

load_data
data = Zobs;

SET.ss   = length(data(1,:)) ;
nobs     = length(data(:,1)) ;
SET.nobs = nobs ;

SET.EST.hidx      = 1:SET.nobs ;
SET.EST.hidx_no_i = 2:SET.nobs ;

SET.EST.obs_ = [
    SET.variable.iobs
    SET.variable.pinfobs
    SET.variable.dy
    SET.variable.dc
    SET.variable.dinve] ;

load ../olg/adj_factors.mat

phi_adj =  phi_s/1 ;
v_adj   =  v_s/1 ;
b_adj   =  A_s/1 ;
a_adj   =  B_s/1 ;
lgv_adj = lgv_vec ;

lgv_adj = [lgv_adj lgv_adj(end)*ones(1,length(year_vec)-length(lgv_adj))];
phi_adj = interp1(year_vec,(1+phi_adj).^1-1,datevec_q,'linear','extrap') ;
v_adj   = interp1(year_vec,(1+v_adj).^1-1,datevec_q,'linear','extrap') ;
b_adj   = interp1(year_vec,(1+b_adj).^1-1,datevec_q,'linear','extrap') ;
a_adj   = interp1(year_vec,(1+a_adj).^1-1,datevec_q,'linear','extrap') ;
tauv    = interp1(year_vec,tauv,datevec_q,'linear','extrap') ;
lgv_adj = interp1(year_vec,(1+lgv_adj).^1-1,datevec_q,'linear','extrap') ;

% restrict changes to a date range
phi_vec = phi_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
v_vec = v_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
b_vec = b_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
a_vec = a_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
tau_vec = tauv(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
lgv_vec = lgv_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;

trend_params = ...
    [phi_vec(1) ; v_vec(1) ; b_vec(1) ; a_vec(1) ; tau_vec(1) ; lgv_vec(1)] ;

save trend_params trend_params

% store vectors for construction of likelihood
SET.tv_str_chg.ant_chg.phi_vec = phi_vec ;
SET.tv_str_chg.ant_chg.v_vec   = v_vec ;
SET.tv_str_chg.ant_chg.b_vec   = b_vec ;
SET.tv_str_chg.ant_chg.a_vec   = a_vec ;
SET.tv_str_chg.ant_chg.tau_vec = tau_vec ;
SET.tv_str_chg.ant_chg.lgv_vec = lgv_vec ;

SET.tv_str_chg.flag = 1 ;

Tbstar = length(SET.tv_str_chg.ant_chg.phi_vec) ;
SET.tv_str_chg.Tbstar = Tbstar ;

SET.tv_str_chg.ant_chg_params   = ...
    {'phi','v','a','b','tau','lgv'} ;
SET.tv_str_chg.unant_chg_params = ...
    {} ; 

gen_set_params(SET,dyn_in) ;

SET.tv_str_chg.cur_t = 1 ;

TT.t_f = zeros(1,length(data)) ;

mats = tvmats(SET,TT) ;

Qt = mats.Qt ; Qt_trendonly = Qt ;
Gt = mats.Gt ; Gt_trendonly = Gt ;

y_ss = (eye(n_-1) - Qt(1:n_-1,1:n_-1,2)) \ Qt(1:n_-1,end,2)  ;
y_ss = [y_ss ; 1] ; % Constant

x_t = y_ss ;
    
for tt=2:Tbstar
    x_t(:,tt) = Qt(:,:,tt) * x_t(:,tt-1) ;%+ Gt(:,:,tt)*n_T(:,tt) ;
end

model_trend = x_t(SET.EST.obs_, 1:end-2) ;

%% Extract Structural Shocks

TT.t_f = T_f ;

mats = tvmats(SET,TT) ;

Qt = mats.Qt ; Qt_trendonly = Qt ;
Gt = mats.Gt ; Gt_trendonly = Gt ;

SET.EST.nparam_est = 11 ; 

SET.tv_str_chg.Qf           = Qt(:,:,SET.trend_dates==(SET.end_date+.25)) ;
SET.tv_str_chg.str_mats     = mats.tv_str_chg.str_mats(:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date)) ;
SET.tv_str_chg.mat_i_f_zlb  = mats.tv_str_chg.mat_i_f_zlb(:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date)) ;
SET.tv_str_chg.zlb          = mats.tv_str_chg.zlb(:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date)) ;
    
Qt = Qt(:,:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date));
Gt = Gt(:,:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date));

V  = zeros(nobs,nobs) ; 
H  = zeros(nobs,n_-1) ; 

Ct = Qt(1:n_-1,end,:) ;
Pt = Qt(1:n_-1,1:n_-1,:) ;
Dt = Gt(1:n_-1,:,:) ;

for ii=1:nobs ; H(ii,SET.EST.obs(ii)) = 1 ; end

[x_T, n_T] = ks_(SET, Ct, Pt, Dt, H, V, data, TT) ;

SET.tv_str_chg.cur_t  = 1 ; 
TT.t_f                = 0*T_f ;

mats = tvmats(SET,TT) ;

Qt_nozlb = mats.Qt ; 
Gt_nozlb = mats.Gt ; 
    
Qt_nozlb = Qt_nozlb(:,:,find(SET.trend_dates==SET.start_date):(SET.maxchk+find(SET.trend_dates==SET.end_date)));
Gt_nozlb = Gt_nozlb(:,:,find(SET.trend_dates==SET.start_date):(SET.maxchk+find(SET.trend_dates==SET.end_date)));

%% Table for Variance Decomposition

table=[] ;

for ii=1:6
    table = [table 100*oo_.conditional_variance_decomposition(:,1,ii)];
end

table = table([...
    SET.variable.iobs
    SET.variable.pinfobs
    SET.variable.dy
    SET.variable.dc
    SET.variable.dinve
    SET.variable.y
    SET.variable.c
    SET.variable.inve],:) ;

disp('') ;
disp('Variance Decomposition, 2Q horizon') ;
latex(table, '%.0f', 'nomath' ) ;

table = [] ;

for ii=1:6
    table = [table 100*oo_.gamma_y{7}(:,ii)];
end

table = table([...
    SET.variable.iobs
    SET.variable.pinfobs
    SET.variable.dy
    SET.variable.dc
    SET.variable.dinve
    SET.variable.y
    SET.variable.c
    SET.variable.inve],:) ;

disp('') ;
disp('Variance Decomposition') ;
latex(table, '%.0f', 'nomath' ) ;

%% IRFs, at 2006, For Appendix

x_t_irf = repmat([(eye(n_-1) - Qt_nozlb(1:n_-1,1:n_-1,find(data_dates==2006))) \ Qt_nozlb(1:n_-1,end,find(data_dates==2006)); 1],[1 100 SET.variable.l_]) ;

for ii=1:SET.variable.l_
    e_t = zeros(SET.variable.l_,100) ;
    e_t(ii,2) = 1 ;
    for tt=2:100
        x_t_irf(:,tt,ii) = Qt_nozlb(:,:,find(data_dates==2006))*x_t_irf(:,tt-1,ii) + ...
            Gt_nozlb(:,:,find(data_dates==2006))*e_t(:,tt) ;
    end
end

figure ;

set(gcf, 'DefaultAxesLineWidth', 1);
set(gcf, 'DefaultLineLineWidth', 1);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',6);
%set(gcf,'PaperPosition',[0 0 11 6]) ;

for ii=1:SET.EST.nobs 
    subplot(3,3,ii) ; hold on ;
    plot(100*(squeeze(x_t_irf(SET.EST.obs(ii),:,:)-1*x_t_irf(SET.EST.obs(ii),1,1)))) ;
    %plot(date_vec,y_non(SET.EST.obs(ii),:),'k','LineWidth',2) ;    
    title(SET.EST.varobs_names{ii},'Interpreter','latex','FontSize',7) ;
    if ii==1
        h=legend('$\eta_\chi$','$\eta_\mu$','$\eta_\xi$','$\eta_R$','$\eta_g$','$\eta_\kappa$','location','northeast');
        set(h,'Interpreter','latex','FontSize',6) ;
    end
    ax = gca ;
    ax.YGrid='on' ;
    xlim([0 12]) ;
end
subplot(3,3,6) ; hold on ;
    plot(1*(squeeze(x_t_irf(SET.variable.y,:,:)-1*x_t_irf(SET.variable.y,1,1)))) ;
    %plot(date_vec,y_non(SET.EST.obs(ii),:),'k','LineWidth',2) ;    
    title('Output','Interpreter','latex','FontSize',7) ;
    ax = gca ;
    ax.YGrid='on' ;
    xlim([0 12]) ;
subplot(3,3,7) ; hold on ;
    plot(1*(squeeze(x_t_irf(SET.variable.c,:,:)-1*x_t_irf(SET.variable.c,1,1)))) ;
    %plot(date_vec,y_non(SET.EST.obs(ii),:),'k','LineWidth',2) ;    
    title('Consumption','Interpreter','latex','FontSize',7) ;
    ax = gca ;
    ax.YGrid='on' ;
    xlim([0 12]) ;
subplot(3,3,8) ; hold on ;
    plot(1*(squeeze(x_t_irf(SET.variable.inve,:,:)-1*x_t_irf(SET.variable.inve,1,1)))) ;
    %plot(date_vec,y_non(SET.EST.obs(ii),:),'k','LineWidth',2) ;    
    title('Investment','Interpreter','latex','FontSize',7) ;
    ax = gca ;
    ax.YGrid='on' ;
    xlim([0 12]) ;
subplot(3,3,9) ; hold on ;
    plot(1*(squeeze(x_t_irf(SET.variable.w,:,:)-1*x_t_irf(SET.variable.w,1,1)))) ;
    %plot(date_vec,y_non(SET.EST.obs(ii),:),'k','LineWidth',2) ;    
    title('Real Wages','Interpreter','latex','FontSize',7) ;
    ax = gca ;
    ax.YGrid='on' ;
    xlim([0 12]) ;

%% Simulations and Counterfactuals

x_t_nodem   = [x_T(:,1) ; 1] ;
x_t_chk     = [x_T(:,1) ; 1] ;
x_t_nozlb   = [x_T(:,1) ; 1] ;
x_t_nosto   = [x_T(:,1) ; 1] ;
x_t_nophish = [x_T(:,1) ; 1] ;
x_t_nopsh   = [x_T(:,1) ; 1] ;
x_t_nomush  = [x_T(:,1) ; 1] ;
x_t_nobsh   = [x_T(:,1) ; 1] ;
x_t_nogsh   = [x_T(:,1) ; 1] ;
x_t_zlb     = [x_T(:,1) ; 1] ;

x_t_nophish_zlb = [x_T(:,1) ; 1] ;
x_t_nopsh_zlb   = [x_T(:,1) ; 1] ;
x_t_nomush_zlb  = [x_T(:,1) ; 1] ;
x_t_nobsh_zlb   = [x_T(:,1) ; 1] ;
x_t_nogsh_zlb   = [x_T(:,1) ; 1] ;

x_t_nophish_nozlb = [x_T(:,1) ; 1] ;
x_t_nopsh_nozlb   = [x_T(:,1) ; 1] ;
x_t_nomush_nozlb  = [x_T(:,1) ; 1] ;
x_t_nobsh_nozlb   = [x_T(:,1) ; 1] ;
x_t_nogsh_nozlb   = [x_T(:,1) ; 1] ;

n_T_nopsh = n_T ;
n_T_nopsh(3,find(data_dates==2008):end) = 0 ;

n_T_nophish = n_T ;
n_T_nophish(1,find(data_dates==2008):end) = 0 ;

n_T_nomush = n_T ;
n_T_nomush(6,find(data_dates==2008):end) = 0 ;

n_T_nobsh = n_T ;
n_T_nobsh(2,find(data_dates==2008):end) = 0 ;

n_T_nogsh = n_T ;
n_T_nogsh(5,find(data_dates==2008):end) = 0 ;

SET.tv_str_chg.flag = 1 ;

SET.maxchk = 40 ;

Qt2 = Qt ;
Gt2 = Gt ;

for tt=2:SET.ss
    
    x_t_nodem(:,tt) = Qt(:,:,1) * x_t_nodem(:,tt-1) + Gt(:,:,1)*n_T(:,tt);
    
    x_t_chk(:,tt) = Qt2(:,:,tt) * x_t_chk(:,tt-1) + Gt2(:,:,tt)*n_T(:,tt);
    
    SET.tv_str_chg.cur_t = tt ;
    
    Qt_in_zlb = Qt_nozlb(:,:,tt:end) ;
    Gt_in_zlb = Gt_nozlb(:,:,tt:end) ;

    SET.tv_str_chg.flag = 1 ;
    
    [x_t_zlb(:,tt), ~, z_f_T(tt)] = ...
        zlb(SET, x_t_zlb(:,tt-1), n_T(:,tt), ...
            Qt_in_zlb, Gt_in_zlb) ;
    
    x_t_nozlb(:,tt) = Qt_nozlb(:,:,tt) * x_t_nozlb(:,tt-1) + ...
        Gt_nozlb(:,:,tt)*n_T(:,tt);    
    
    x_t_nosto(:,tt) = Qt_nozlb(:,:,tt) * x_t_nosto(:,tt-1) ;
    
    % For counterfactuals without shocks
    if data_dates(tt)>2008.75
        
    SET.tv_str_chg.flag = 0 ;
    Qt_in_zlb(:,:,1:SET.maxchk) = ...
        repmat(Qt_trendonly(:,:,find(SET.trend_dates==2008.75)),[1 1 SET.maxchk]) ;
    Gt_in_zlb(:,:,1:SET.maxchk) = ...
        repmat(Gt_trendonly(:,:,find(SET.trend_dates==2008.75)),[1 1 SET.maxchk]) ;

    SET.zlb.Qf          = Qt_trendonly(:,:,find(SET.trend_dates==2008.75)) ; 
    SET.zlb.mat_init    = mats.tv_str_chg.str_mats(:,find(SET.trend_dates==2008.75)) ; SET.zlb.mat_init = SET.zlb.mat_init{1} ;
    SET.zlb.mat_i_f_zlb = mats.tv_str_chg.mat_i_f_zlb(:,find(SET.trend_dates==2008.75)) ; SET.zlb.mat_i_f_zlb = SET.zlb.mat_i_f_zlb{1} ;
    
    end
    
    [x_t_nopsh(:,tt), ~, T_Z(tt,1)] = ...
        zlb(SET, x_t_nopsh(:,tt-1), n_T_nopsh(:,tt), ...
            Qt_in_zlb, Gt_in_zlb) ;
        
    [x_t_nobsh(:,tt), ~, T_Z(tt,2)] = ...
        zlb(SET, x_t_nobsh(:,tt-1), n_T_nobsh(:,tt), ...
            Qt_in_zlb, Gt_in_zlb) ;        
        
    [x_t_nomush(:,tt), ~, T_Z(tt,3)] = ...
        zlb(SET, x_t_nomush(:,tt-1), n_T_nomush(:,tt), ...
            Qt_in_zlb, Gt_in_zlb) ;        
        
    [x_t_nophish(:,tt), ~, T_Z(tt,4)] = ...
        zlb(SET, x_t_nophish(:,tt-1), n_T_nophish(:,tt), ...
            Qt_in_zlb, Gt_in_zlb) ;        
        
    [x_t_nogsh(:,tt), ~, T_Z(tt,5)] = ...
        zlb(SET, x_t_nogsh(:,tt-1), n_T_nogsh(:,tt), ...
            Qt_in_zlb, Gt_in_zlb) ;        
    
    x_t_nopsh_zlb(:,tt) = Qt2(:,:,tt) * x_t_nopsh_zlb(:,tt-1) + ...
        Gt2(:,:,tt) * n_T_nopsh(:,tt) ;
    
    x_t_nobsh_zlb(:,tt) = Qt2(:,:,tt) * x_t_nobsh_zlb(:,tt-1) + ...
        Gt2(:,:,tt) * n_T_nobsh(:,tt) ;
    
    x_t_nomush_zlb(:,tt) = Qt2(:,:,tt) * x_t_nomush_zlb(:,tt-1) + ...
        Gt2(:,:,tt) * n_T_nomush(:,tt) ;
    
    x_t_nophish_zlb(:,tt) = Qt2(:,:,tt) * x_t_nophish_zlb(:,tt-1) + ...
        Gt2(:,:,tt) * n_T_nophish(:,tt) ;
    
    x_t_nogsh_zlb(:,tt) = Qt2(:,:,tt) * x_t_nogsh_zlb(:,tt-1) + ...
        Gt2(:,:,tt) * n_T_nogsh(:,tt) ;

	x_t_nopsh_nozlb(:,tt) = Qt_nozlb(:,:,tt) * x_t_nopsh_nozlb(:,tt-1) + ...
        Gt_nozlb(:,:,tt) * n_T_nopsh(:,tt) ;
    
    x_t_nobsh_nozlb(:,tt) = Qt_nozlb(:,:,tt) * x_t_nobsh_nozlb(:,tt-1) + ...
        Gt_nozlb(:,:,tt) * n_T_nobsh(:,tt) ;
    
    x_t_nomush_nozlb(:,tt) = Qt_nozlb(:,:,tt) * x_t_nomush_nozlb(:,tt-1) + ...
        Gt_nozlb(:,:,tt) * n_T_nomush(:,tt) ;
    
    x_t_nophish_nozlb(:,tt) = Qt_nozlb(:,:,tt) * x_t_nophish_nozlb(:,tt-1) + ...
        Gt_nozlb(:,:,tt) * n_T_nophish(:,tt) ;
    
    x_t_nogsh_nozlb(:,tt) = Qt_nozlb(:,:,tt) * x_t_nogsh_nozlb(:,tt-1) + ...
        Gt_nozlb(:,:,tt) * n_T_nogsh(:,tt) ;
    
end

%% Figure 6: Keep Demographics Fixed at Certain Dates

sh_irf = zeros(SET.variable.l_,40) ;
sh_irf(SET.variable.shock.eps_mu,2) = 10 ; % Large financial shock

Qt_in_1990 = repmat(Qt_trendonly(:,:,find(SET.trend_dates==1990)),[1 1 200]) ;
Gt_in_1990 = repmat(Gt_trendonly(:,:,find(SET.trend_dates==1990)),[1 1 200]) ;

Qt_in_2008 = repmat(Qt_trendonly(:,:,find(SET.trend_dates==2008)),[1 1 200]) ;
Gt_in_2008 = repmat(Gt_trendonly(:,:,find(SET.trend_dates==2008)),[1 1 200]) ;

x_t_irf_1990 = [(eye(n_-1) - Qt_in_1990(1:n_-1,1:n_-1,1)) \ Qt_in_1990(1:n_-1,end,1) ; 1] ;
x_t_irf_2008 = [(eye(n_-1) - Qt_in_2008(1:n_-1,1:n_-1,1)) \ Qt_in_2008(1:n_-1,end,1) ; 1] ;

SET.tv_str_chg.flag = 0 ;
SET.zlb.Qf          = Qt_trendonly(:,:,find(SET.trend_dates==1990)) ; 
SET.zlb.mat_init    = mats.tv_str_chg.str_mats(:,find(SET.trend_dates==1990)) ; SET.zlb.mat_init = SET.zlb.mat_init{1} ;
SET.zlb.mat_i_f_zlb = mats.tv_str_chg.mat_i_f_zlb(:,find(SET.trend_dates==1990)) ; SET.zlb.mat_i_f_zlb = SET.zlb.mat_i_f_zlb{1} ;

for tt=2:40
    
[x_t_irf_1990(:,tt), ~, ~] = ...
        zlb(SET, x_t_irf_1990(:,tt-1), sh_irf(:,tt), ...
            Qt_in_1990, Gt_in_1990) ;
        
end

SET.tv_str_chg.flag = 0 ;
    SET.zlb.Qf          = Qt_trendonly(:,:,find(SET.trend_dates==2008)) ; 
    SET.zlb.mat_init    = mats.tv_str_chg.str_mats(:,find(SET.trend_dates==2008)) ; SET.zlb.mat_init = SET.zlb.mat_init{1} ;
    SET.zlb.mat_i_f_zlb = mats.tv_str_chg.mat_i_f_zlb(:,find(SET.trend_dates==2008)) ; SET.zlb.mat_i_f_zlb = SET.zlb.mat_i_f_zlb{1} ;

for tt=2:40
        
[x_t_irf_2008(:,tt), ~, ~] = ...
        zlb(SET, x_t_irf_2008(:,tt-1), sh_irf(:,tt), ...
            Qt_in_2008, Gt_in_2008) ;        
        
end

figure; 

set(gcf, 'DefaultAxesLineWidth', 1);
set(gcf, 'DefaultLineLineWidth', 1.5);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',7.5);
%set(gcf,'PaperPosition',[0 0 11 6]) ;

subplot(2,2,1) ; hold on ;
    plot(400*x_t_irf_1990(SET.variable.iobs,:))
    plot(400*x_t_irf_2008(SET.variable.iobs,:),'-.')
    xlim([0 20]) ;
    
    title('A. Fed Funds Rate','Interpreter','latex') ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
    
    h = legend('Under 1990 Demographics','Under 2008 Demographics', 'Location','NorthEast');
    set(h,'Interpreter','latex','FontSize',7.5);      
    
subplot(2,2,2) ; hold on ;
    plot(100*log(x_t_irf_1990(SET.variable.y,:)./x_t_irf_1990(SET.variable.y,1)))
    plot(100*log(x_t_irf_2008(SET.variable.y,:)./x_t_irf_2008(SET.variable.y,1)),'-.')
    
    line([0 20],[0 0],'Color','k','LineWidth',1) ;
    xlim([0 20]) ;
    ylim([-8 2]) ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;

    title('B. Output','Interpreter','latex') ; 

%% Figure 1: Output, Interest Rates, and the Employment-Pop Ratio

data_for_figs

gr_rate_86_to_2008 = mean(x_T(SET.variable.dy,find(data_dates==1986):...
    find(data_dates==2007.75))) ;

ypredict      = log(gen_idx(x_t(SET.variable.dy,:)*0+gr_rate_86_to_2008)) ;  
yplotdat      = log(gen_idx(x_T(SET.variable.dy,:))) ;  
yplot_1984dem = log(gen_idx(x_t_nodem(SET.variable.dy,:))) ;  
yplot_zlb     = log(gen_idx(x_t_zlb(SET.variable.dy,:))) ;  
yplot         = log(gen_idx(x_t(SET.variable.dy,:))) ;  
yplot_nozlb   = log(gen_idx(x_t_nozlb(SET.variable.dy,:))) ;  

base = 2007 ;

figure;

set(gcf, 'DefaultAxesLineWidth', .75);
set(gcf, 'DefaultLineLineWidth', 1.25);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',7);
set(gcf, 'Position', [1.6353e+03 712 560*(2) 420*(2/3)])

subplot(1,3,1) ; hold on;
    
    plot(data_dates,yplotdat-yplotdat(find(data_dates==base))+1) ;
    plot(trend_dates,ypredict-ypredict(find(trend_dates==base))+1,'--k','LineWidth',.75);  
    
    xlim([1990 2021]) ;
        
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
	
    h = legend('','Trend, 1986 to 2007', 'Location','northwest');
    set(h,'Interpreter','latex','FontSize',6);    
    
    title('A. Output (2007 = 1)','Interpreter','latex') ;
    
subplot(1,3,2) ; hold on;
    plot(data_dates,100*(exp(x_T(SET.variable.iobs,:)).^4-1));  
    plot(rir_treasury_q_dates,rir_treasury_q,'-.') ; % 10 Year Treasury minus Core Inflation
    xlim([1990 2021]) ;
	line([1986 2030],[0 0],'LineWidth',0.5,'Color','k') ;

    ylim([-2 10]) ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;

    title('B. Interest Rates, Annual \%','Interpreter','latex') ;
    
    h = legend('Fed Funds Rate','10Y Treasury - Core Infl.','', 'Location','northeast');
    set(h,'Interpreter','latex','FontSize',5);  
    
subplot(1,3,3) ; hold on ;
    plot(emp_pop_ratio_dates,emp_pop_ratio)
    xlim([1960 2060]) ;
    title('C. Employment to Population, \%','Interpreter','latex')
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
    xlim([1990 2021]) ;    
    ylim([50 70]) ;   

%% Figure 3, Demographic trends only

figure; 

set(gcf, 'DefaultAxesLineWidth', .75);
set(gcf, 'DefaultLineLineWidth', 1.25);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',7);
set(gcf, 'Position', [1.6353e+03 712 560*(2) 420*(2/3)])

subplot(1,3,2) ; hold on;
    plot(data_dates,100*(exp(x_T(SET.variable.iobs,:)).^4-1));  
	plot(trend_dates,100*(exp(x_t(SET.variable.iobs,:)).^4-1),'-.');  
    xlim([1990 2030]) ;

    ylim([0 10]) ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;

    title('B. Fed Funds Rate, Annual \%','Interpreter','latex') ;
    
subplot(1,3,1) ; hold on;
    
    plot(data_dates,yplotdat-yplotdat(find(data_dates==base))+1) ;
	plot(trend_dates,yplot-yplot(find(trend_dates==base))+1,'-.') ;
    plot(trend_dates,ypredict-ypredict(find(trend_dates==base))+1,'--k','LineWidth',1);  
    
    xlim([1990 2030]) ;
        
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
	
    h = legend('','','Trend, 1986 to 2007', 'Location','northwest');
    set(h,'Interpreter','latex','FontSize',6);    
    
    title('A. Output (2007 = 1)','Interpreter','latex') ;
subplot(1,3,3) ; hold on;
    plot(data_dates,100*((1+x_T(SET.variable.rir,:)).^4-1));  
    plot(trend_dates,100*((1+x_t(SET.variable.rir,:)).^4-1),'-.');  

    line([1986 2030],[0 0],'LineWidth',0.5,'Color','k') ;
    xlim([1990 2030]) ;
    
    h = legend('$i_t-E_t\pi_{t+1}$', 'Location','northeast');
    set(h,'Interpreter','latex','FontSize',7);  
    
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
    
    title('C. Real Interest Rate, Annual \%','Interpreter','latex') ;

%% Figure 5, Demographics Fixed from 1986

figure; 

set(gcf, 'DefaultAxesLineWidth', .75);
set(gcf, 'DefaultLineLineWidth', 1.25);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',7);
set(gcf, 'Position', [1.6353e+03 712 560*(2) 420*(2/3)])

subplot(1,3,2) ; hold on;
    plot(data_dates,100*(exp(x_T(SET.variable.iobs,:)).^4-1));  
	plot(data_dates,100*(exp(x_t_nodem(SET.variable.iobs,:)).^4-1),'-.');  
    xlim([1990 2020]) ;

    ylim([0 10]) ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;

    title('B. Fed Funds Rate, Annual \%','Interpreter','latex') ;
    
subplot(1,3,1) ; hold on;
    
    plot(data_dates,yplotdat-yplotdat(find(data_dates==base))+1) ;
	plot(data_dates,yplot_1984dem-yplot_1984dem(find(data_dates==base))+1,'-.') ;
    plot(trend_dates,ypredict-ypredict(find(trend_dates==base))+1,'--k','LineWidth',1);  
    
    xlim([1990 2020]) ;
        
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
	
    h = legend('','','Trend, 1986 to 2007', 'Location','northwest');
    set(h,'Interpreter','latex','FontSize',6);    
    
    title('A. Output (2007 = 1)','Interpreter','latex') ;
subplot(1,3,3) ; hold on;
    plot(data_dates,100*((1+x_T(SET.variable.rir,:)).^4-1));  
    plot(data_dates,100*((1+x_t_nodem(SET.variable.rir,:)).^4-1),'-.');  
    
    line([1986 2030],[0 0],'LineWidth',0.5,'Color','k') ;
    xlim([1990 2020]) ;
    
    h = legend('$i_t-E_t\pi_{t+1}$', 'Location','northeast');
    set(h,'Interpreter','latex','FontSize',6);  
    
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
    
    title('C. Real Interest Rate, Annual \%','Interpreter','latex') ;

%% Figure 7: No Forward Guidance

figure; 

set(gcf, 'DefaultAxesLineWidth', .75);
set(gcf, 'DefaultLineLineWidth', 1.25);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',7);
set(gcf, 'Position', [1.6353e+03 712 560*(2) 420*(2/3)])

subplot(1,3,2) ; hold on;
    plot(data_dates,100*(exp(x_T(SET.variable.iobs,:)).^4-1));  
	plot(data_dates,100*(exp(x_t_zlb(SET.variable.iobs,:)).^4-1),'-.');  

    xlim([1990 2020]) ;

    ylim([0 10]) ;
    ax = gca ;
    ax.YGrid='on';
    box('off') ;

    title('B. Fed Funds Rate, Annual \%','Interpreter','latex') ;
    
subplot(1,3,1) ; hold on;
    
    plot(data_dates,yplotdat-yplotdat(find(data_dates==base))+1) ;
	plot(data_dates,yplot_zlb-yplot_zlb(find(data_dates==base))+1,'-.') ;
    plot(trend_dates,ypredict-ypredict(find(trend_dates==base))+1,'--k','LineWidth',1);  

    xlim([1990 2020]) ;
        
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
	
    h = legend('','','Trend, 1986 to 2007', 'Location','northwest');
    set(h,'Interpreter','latex','FontSize',6);    
    
    title('A. Output (2007 = 1)','Interpreter','latex') ;
subplot(1,3,3) ; hold on;
    plot(data_dates,100*((1+x_T(SET.variable.rir,:)).^4-1));  
    plot(data_dates,100*((1+x_t_zlb(SET.variable.rir,:)).^4-1),'-.');  

    line([1986 2030],[0 0],'LineWidth',0.5,'Color','k') ;
    xlim([1990 2020]) ;
    
    h = legend('$i_t-E_t\pi_{t+1}$', 'Location','northeast');
    set(h,'Interpreter','latex','FontSize',6);  
    
    ax = gca ;
    ax.YGrid='on';
    box('off') ;
    
    title('C. Real Interest Rate, Annual \%','Interpreter','latex') ;

%% Figure 8

gr_rate_86_to_2008 = mean(x_T(SET.variable.dy,find(data_dates==1986):...
    find(data_dates==2007.75))) ;

date = date_vec;

x_t_nosto_nozlb = x_t_nosto ;

y_idx       = gen_idx(x_T(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nozlb = gen_idx(x_t_nozlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_zlb   = gen_idx(x_t_zlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_f     = gen_idx(x_T(SET.variable.dy_f,:),1:1:length(x_T(1,:))) ;
y_idx_nosto = gen_idx(x_t_nosto_nozlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_trend = gen_idx(0*x_t_nosto_nozlb(SET.variable.dy,:)+gr_rate_86_to_2008,1:1:length(x_T(1,:))) ;
y_idx_nodem = gen_idx(x_t(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nopsh = gen_idx(x_t_nopsh(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nophish = gen_idx(x_t_nophish(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nobsh = gen_idx(x_t_nobsh(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nogsh = gen_idx(x_t_nogsh(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nomush = gen_idx(x_t_nomush(SET.variable.dy,:),1:1:length(x_T(1,:))) ;

y_idx_nopsh_zlb = gen_idx(x_t_nopsh_zlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nophish_zlb = gen_idx(x_t_nophish_zlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nobsh_zlb = gen_idx(x_t_nobsh_zlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nogsh_zlb = gen_idx(x_t_nogsh_zlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nomush_zlb = gen_idx(x_t_nomush_zlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;

y_idx_nopsh_nozlb = gen_idx(x_t_nopsh_nozlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nophish_nozlb = gen_idx(x_t_nophish_nozlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nobsh_nozlb = gen_idx(x_t_nobsh_nozlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nogsh_nozlb = gen_idx(x_t_nogsh_nozlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;
y_idx_nomush_nozlb = gen_idx(x_t_nomush_nozlb(SET.variable.dy,:),1:1:length(x_T(1,:))) ;

base = 2007 ;

% Panel A: Decomposition showing effect of demographics
figure; 

set(gcf, 'DefaultAxesLineWidth', 1);
set(gcf, 'DefaultLineLineWidth', 2.5);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',14);

hold on ;

plot(date_vec,100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx./y_idx(find(date==base)))))

plot(date_vec,100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx_zlb./y_idx_zlb(find(date==base)))),'-.')

plot(date_vec,100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx_nosto./y_idx_nosto(find(date==base)))),'--')

xlim([2007 2020]) ;
line([2000 2020],[0 0],'Color','k','LineWidth',1) ;
ylim([-20 5]) ;

ax = gca ;
ax.YGrid='on';
box('off') ;

h=legend('Data and Model','No Forward Guidance','Demographics Only','location','southwest') ;
set(h,'Interpreter','latex','FontSize',8.5);  

% Panel B: Shock decomposition without ZLB
figure; 

set(gcf, 'DefaultAxesLineWidth', 1);
set(gcf, 'DefaultLineLineWidth', 2.5);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',14);

hold on ;

plot(date_vec,100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx./y_idx(find(date==base)))))

tmp = 100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx./y_idx(find(date==base)))) ;

tmp2 = 100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx_nomush_nozlb./y_idx_nomush_nozlb(find(date==base)))) ;

 plot(date_vec,tmp-tmp2,'-x','MarkerSize', 9, 'MarkerIndices',1:4:length(date_vec))

 tmp2 = 100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx_nopsh_nozlb./y_idx_nopsh_nozlb(find(date==base)))) ;

 plot(date_vec,tmp-tmp2,'-o','MarkerSize', 9, 'MarkerIndices',1:4:length(date_vec))

 tmp2 = 100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx_nobsh_nozlb./y_idx_nobsh_nozlb(find(date==base)))) ;

 plot(date_vec,tmp-tmp2,'--')

 tmp2 = 100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx_nogsh_nozlb./y_idx_nogsh_nozlb(find(date==base)))) ;

 plot(date_vec,tmp-tmp2,'-.')     

tmp2 = 100*(-log(y_idx_trend./y_idx_trend(find(date==base)))...
    +log(y_idx_nophish_nozlb./y_idx_nophish_nozlb(find(date==base)))) ;

 plot(date_vec,tmp-tmp2,'-.*','MarkerSize', 9, 'MarkerIndices',1:4:length(date_vec))       

xlim([2007 2020]) ;
line([2000 2020],[0 0],'Color','k','LineWidth',1) ;
ylim([-20 5]) ;

ax = gca ;
ax.YGrid='on';
box('off') ;

h=legend('Data and Model','Investment','Markup','Technology','Government','Preference','location','southwest') ;
set(h,'Interpreter','latex','FontSize',10);  

%% Statistics in Paper

tmp1=100*(-log(y_idx_trend./y_idx_trend(find(date==base)))+log(y_idx./y_idx(find(date==base)))) ;
    
disp('') ;
disp('Gap Between US Output PC and Trend, 2019Q4:') ;
disp(tmp1(find(date_vec==2019.75)))
disp('') ;

tmp2 = 100*(-log(y_idx_trend./y_idx_trend(find(date==base)))+log(y_idx_nosto./y_idx_nosto(find(date==base)))) ;

disp('') ;
disp('Gap Between US Output PC and Trend, Demographics 2019Q4:') ;
disp(tmp2(find(date_vec==2019.75)))
disp('') ;
    